A multi-omic dissection of super-enhancer driven oncogenic gene expression programs in ovarian cancer

The human genome contains regulatory elements, such as enhancers, that are often rewired by cancer cells for the activation of genes that promote tumorigenesis and resistance to therapy. This is especially true for cancers that have little or no known driver mutations within protein coding genes, such as ovarian cancer. Herein, we utilize an integrated set of genomic and epigenomic datasets to identify clinically relevant super-enhancers that are preferentially amplified in ovarian cancer patients. We systematically probe the top 86 super-enhancers, using CRISPR-interference and CRISPR-deletion assays coupled to RNA-sequencing, to nominate two salient super-enhancers that drive proliferation and migration of cancer cells. Utilizing Hi-C, we construct chromatin interaction maps that enable the annotation of direct target genes for these super-enhancers and confirm their activity specifically within the cancer cell compartment of human tumors using single-cell genomics data. Together, our multi-omic approach examines a number of fundamental questions about how regulatory information encoded into super-enhancers drives gene expression networks that underlie the biology of ovarian cancer.

O varian cancer is one of the deadliest cancers among women worldwide and is the leading cause of gynecologic-related cancer deaths in the US 1 . High-grade serous ovarian cancer (HGSOC) is the most common subtype (~80% of all ovarian cancer) and is characterized by a high number of copy number alterations and few driver mutations, which is thought to account for the clinical aggressiveness of this disease as well as the eventual development of chemoresistance 2,3 . The most commonly seen mutation in HGSOC is p53 (>90% of cases), followed by a low, but statistically significant, prevalence of recurrent somatic mutations in NF-1, BRCA1/2, and CDK2, which often lead to genomic instability [4][5][6] . Due to this genomic instability, ovarian cancer has a high rate of copy number abnormalities and recent studies have shown that these alterations can be used to stratify HGSOC 2 . However, the paucity of known driver mutations for ovarian cancer has made it difficult to develop effective targeted therapies. Consequently, the standard of care remains cytoreductive surgery followed by carboplatin/taxane chemotherapy, with~75% of patients experiencing a recurrence 2,3 . Thus, additional analysis of the non-coding regions of the genome, that extends beyond gene profiling, is desperately needed.
Mounting evidence suggests that regulatory elements, such as transcriptional enhancers, can be rewired or hijacked by cancer cells for the activation of genes that promote tumor formation, metastasis, and resistance to therapy [7][8][9] . This is especially true for cancers that have little or no known driver mutations within protein-coding genes, such as ovarian cancer 10 . Enhancers are non-coding DNA elements that contain information for the binding of transcription factors and interact spatially with their target genes to orchestrate spatiotemporal patterns of gene expression 11,12 . It is estimated that there are hundreds of thousands of enhancers found throughout our genome and these can act independent of orientation and linear distance from their target genes, forming high-order chromatin loops with their target genes. Of note, the activity of enhancers is often restricted to a particular cell type or to specific physiological or pathological conditions, enabling their genomic function to determine precisely when, where, and at what level each of our genes is expressed [13][14][15] . Large clusters of neighboring enhancers that have unusually high occupancy of interacting factors are typically called super-enhancers (SEs) 16 . These super-enhancers are known to regulate key cell identity genes, and in cancer are known to drive oncogene expression 17 .
The high transcriptional output of cancer cells is thought to be sustained by the activity of super-enhancers, suggesting cancer cells can become addicted to super-enhancer-driven regulatory networks 18 . Furthermore, recent studies in ovarian cancer have demonstrated the capacity of super-enhancers and their associated networks of transcription factors to directly influence chemoresistance 19,20 . The molecular characteristics and high activity of super-enhancers make them exquisitely sensitive to epigenetic drugs, more so than typical enhancers 21 . Thus, there is a growing belief that exploiting transcriptional dependence by targeting oncogenic super-enhancers may be a valid therapeutic avenue 21 . For example, the bromodomain-containing protein 4 (BRD4) is a druggable transcription factor that recognizes acetylated histone proteins and is found in large quantities at super-enhancers 22,23 . Small molecule inhibition of BRD4 (such as JQ1 and BET inhibitors) has been shown to reduce cell proliferation and survival in vivo as well as increase therapeutic sensitivity of several cancer types, leading to the development of several clinical trials 20,23,24 . However, despite their effectiveness in inhibiting oncogenic processes in ovarian cancer cells, anti-BRD4 agonists remain a poor therapeutic option due to their overall toxicity and delivery constraints 25 . Nevertheless, the study of BRD4-associated super-enhancers in ovarian cancer may lead to the identification of biomarkers, downstream druggable targets, and a better understanding of the regulatory processes that drive this disease.
To this end, the studies described herein examine several fundamental questions about how regulatory information is encoded into super-enhancers, how they are preferentially amplified in ovarian cancer cells, and how they drive gene expression networks that underlie the biology of ovarian cancer cells. We use an integrated genomic and computational framework to (1) identify BRD4-enriched and copy number amplified super-enhancers in ovarian cancer patients, (2) systematically probe the functions of the top 86 ovarian cancer-enriched superenhancers using CRISPR interference assays (CRISPRi) (dCas9-KRAB) coupled to RNA-seq, (3) validate their roles in driving the proliferation and migration of cancer cells via CRISPR-knockouts, (4) annotate direct target genes using chromatin looping information via Hi-C, and (5) confirm their activity specifically within the cancer cell compartment of human tumors using single-cell genomics data.

Results
Identification of BRD4-enriched super-enhancers in ovarian cancer. Super-enhancers are one of the most salient regulatory elements in the genome and are known to be repurposed by cancer cells to drive the expression of oncogenes 16,26 . Due to the unusually high levels of interacting transcription factors and the prominence of their target genes, super-enhancers contain untapped potential that can lead to a new set of markers with diagnostic and prognostic potential, or even serve as tractable targets for therapeutic intervention 21,27 . To identify enhancers likely to be associated with oncogenic gene expression programs, we leveraged both ovarian cancer cell line epigenetic data and patient tumor RNA-seq and copy number data from The Cancer Genome Atlas (TCGA) 10 (Fig. 1a).
First, we used existing ChIP-seq data in the well-vetted highgrade serous ovarian cancer cell line OVCAR3 to identify active enhancers by searching for co-localization of the histone modification histone H3 lysine 27 acetylation (H3K27ac) and BRD4 (Fig. 1f) [28][29][30] . BRD4 enrichment was considered a critical component for the detection of potentially oncogenic enhancers due to key observations previously shown in ovarian cancer patients 23 . Namely, across the entirety of the TCGA Pan-Cancer dataset, ovarian cancer patients have the highest rate of genetic amplifications at the BRD4 locus, with~11% of patients having an amplification of this region (Fig. 1b) [5][6][7]31 . Moreover, ovarian cancer has the highest overall expression of BRD4 across all TCGA cancer types and patients with increased expression of BRD4 experienced significantly reduced survival times as determined through Kaplan-Meier analysis (Fig. 1c, d) 6,[31][32][33][34][35] . Therefore, we defined active enhancers as intergenic regions that contained at least a 1-base pair overlap between statistically significant BRD4 peaks and H3K27ac peaks called by the MACS2 peak calling algorithm (Fig. 1e) 36 . To focus on distal enhancer elements, any peaks that overlapped with annotated genes or promoter regions were removed. This pipeline identified 12,339 BRD4-enriched active enhancer elements in ovarian cancer cells.
To determine if these enhancers are lineage-specific or extensible to other cancer types, we investigated the overlap with existing enhancer annotations across all normal tissues (defined by the ENCODE consortium), and across existing annotations in other cancer types ( Supplementary Fig. 1c, d) 16,37 . In addition to this general analysis, we also specifically compared our 12,339 enhancer predictions to previously annotated enhancers in two distinct fallopian tube secretory epithelial cell lines (FTSEC) 38 , currently thought to be the most common precursor cell of origin for HGSOC ( Supplementary Fig. 1b) 39,40 . We found that 44.1% of the 12,339 BRD4-enriched enhancers had at least 1-base pair overlap with active enhancers in normal tissues and this number increases to 73.6% when compared to active enhancers across several cancer types. Of particular interest, we also identified 6000 enhancers in HGSOC cells that are not active in normal FTSEC cells (Supplementary Fig. 1b). The aforementioned importance of BRD4 and the high degree of overlap between the enhancers identified in this study with previously annotated enhancers in cancer cells gave us confidence for using these data for calling super-enhancers.
From our pool of 12,339 constituent enhancers, we identified 126 super-enhancer regions using the rank ordering of superenhancers (ROSE) algorithm (Fig. 1g, h, and Supplementary Data 1-4) 41,42 . To determine if these BRD4-enriched superenhancers are relevant to ovarian cancer patients, we leveraged the assay for transposase-accessible chromatin at single-cell resolution (scATAC-seq) data generated from HGSOC patients to measure the activity of the super-enhancers within these   Fig. 1a). Taken together, these data suggest that the super-enhancers identified using our pipeline are not cell line specific and may be relevant to both ovarian cancer and other cancer types. To further investigate the clinical utility of these SEs, we next looked for evidence in patient tumors using both TCGA RNA-seq and copy number variation data.

Copy Number Variation and expression Quantitative Trait
Loci (CNVeQTL) analysis nominates potentially oncogenic super-enhancers. Given that copy number variation (CNV) has been previously identified as an important hallmark of ovarian cancer, we sought to investigate whether these BRD4-enriched super-enhancers were preferentially amplified in ovarian cancer patients 2 . To this end, we performed a computational experiment making use of publicly available copy number variation data across~600 ovarian cancer patients 10 to compare the copy number amplification values overlapping our SE regions to the amplification across the ovarian cancer genome as a whole, by both random-draw (pseudo-bootstrap) and direct comparison analyses (Fig. 2a). Copy number variation values across~600 ovarian cancer patients were quantified by dividing the genome into uniform 15 kb sliding windows and assigning CNV segment values within each window (Fig. 2d). We then compared the amplification of the windows that overlap the SEs against an equivalent number of randomly drawn windows across the ovarian cancer genome (inclusive of our SE regions). The random drawing of windows was iterated 10,000 times and, in each comparison, there was significant enrichment in amplification of the SE overlapping windows compared to the random groups (Fig. 2c). This observation was reinforced by comparing SE CNV to the CNV across the ovarian cancer genome as a whole (Fig. 2e). Remarkably, amplification of the super-enhancers themselves was prognostic of clinical outcome 44 . In many cases, patients with increased copy numbers had significantly increased hazard ratio and reduced survival times, suggesting that superenhancer copy numbers may be of prognostic value (Fig. 2b). Taken together, these data suggest that the SEs we identified in OVCAR3 cells are preferentially amplified in ovarian cancer patients and that some SE amplifications are associated with reduced survival.
To better understand how amplification of these SEs are associated with oncogenic gene expression networks, we leveraged the RNA-seq data generated from a subset of the same ovarian cancer patients (~300) to link the SEs to gene expression. We took inspiration from a commonly used approach in complex genetics which associates nucleotide variants to changes in gene expression called eQTL analysis 45,46 . However, unlike eQTL analysis which focuses on point mutations, the comparison, in this case, focuses on changes in copy number across SE loci to changes in gene expression within each patient (copy number variation expression quantitative trait loci (CNVeQTL)) (Supplementary Fig. 2a). The assumption is that amplification or deletion of SE regions should affect their target genes, therefore, looking across hundreds of patients for shared patterns of variation will identify putative target genes of each SE. However, since we altered the input data of the eQTL detection software to utilize two quantitative variables (copy number and gene expression), we needed to determine a robust indication of our null condition for statistical analysis.
To generate the null dataset, we broke the linkage of RNA to copy number by randomly permuting the columns of the RNA data matrix and then running Matrix eQTL 46 on this permutated dataset, repeating this process 100k times, and using the median distribution from all 100k trials to inform our experimental analysis. Importantly, all 100k runs using the permutated null data showed a relatively uniform distribution of p-values across the null condition, suggesting no meaningful relationship between copy number and gene expression, and returned a similar count of total significant CNVeQTLs (the median number of CNVeQTLs across all 100k was 11,632) ( Supplementary  Fig. 2b). In contrast, the results from the true data show a much sharper peak around p-value = 0 and returned a much larger number of significant CNVeQTLs (n = 126,438) (Supplementary Fig. 2c and Supplementary Data 5). We used the results of the 100k null experiments to determine an empirical false discovery rate 47 of about 0.092. This data also allowed us to investigate some higher-order questions, such as whether the number of CNVeQTL detected was strictly a function of size. While there was a modest linear relationship between these features, this analysis suggested something other than genomic size influenced the number of CNVeQTL ( Supplementary Fig. 2d). Collectively, these data suggest that amplification of the super-enhancer regions is associated with pervasive gene expression changes in human tumors, reinforcing the idea they are not merely cell line specific, and they may be preferentially amplified for a biologically meaningful reason.
We recognize that the identification of 126,438 CNVeQTL linkages across 126 super-enhancers seems high, despite the null distributions tested, and that the vast majority of copy number amplifications will have very strong effects in cis (and most will have effects in trans) irrespective of their designation as a super-enhancer. Therefore, to functionally validate and assess the full scope of this data, we chose the top 86 super-enhancers ranked by BRD4 enrichment and H3K27ac signal (which were located both above and below the CNVeQTL prediction line) to perturb using a high throughput CRISPRi screen ( Supplementary  Fig. 2d). Fig. 1 Identification of BRD4-enriched super-enhancers in ovarian cancer. a Flowchart of the analysis strategy used to identify clinically relevant BRD4enriched SEs in ovarian cancer. b Bar chart depicting the alteration frequency of the BRD4 locus across the top 16 highest altered cancer types in the TGCA Pan-Cancer patient cohort (ovarian cancer = OV) retrieved from cBioPortal 6,31 . c Box plots showing normalized BRD4 expression across the top 16 highest expressing cancer types in the TCGA Pan-Cancer patient cohort (ovarian cancer = OV) retrieved from cBioPortal 6,31 . Boxplot is centered on the median (center line), with the upper and lower quartiles creating the bounds of the box (the IQR). The minimum and maximum values, after disregarding outliers, are represented by the upper and lower whiskers. d Kaplan-Meier plots 34,35 showing the relationship between BRD4 expression and progression-free survival in ovarian cancer patients with high-grade serous (n = 1232) or endometrioid histology (n = 62). Patients are split by median expression of BRD4. The red line represents patients in the high expression cohort and the black line low expression cohort. e Cartoon depicting the analysis strategy for integrating H3K27ac and BRD4 ChIP-seq data and selecting overlapping peaks to call super-enhancers. BRD4 is shown in green and H3K27ac in blue. f Top: Meta-ChIP plot of the signal across shared peaks showing overlap of H3K27ac and BRD4 signal. Bottom: Heatmap of ChIP signal across all 12,339 called shared peaks. The samples are scaled relative to the background for that signal group independent of the other signal (BRD4 to BRD4 background; H3K27ac to H3K27ac background). g BRD4 signal versus enhancer rank plot showing the identification of 126 super-enhancers as defined by the ROSE software. h Tabulation of the total number of enhancers/peaks identified. Source data are provided as a Source Data file.
High throughput CRISPR-interference screen highlights super-enhancer target gene relationships. To systematically probe the functions of each SE and determine the consequences on gene expression, we used high-throughput CRISPR-interference assays coupled to RNA-seq. For this experiment, we engineered OVCAR3 cells to stably express nuclease deficient Cas9 fused to the KRAB effector domain (dCas9-KRAB). The KRAB effector domain induces local chromatin repression via methylation of histone 3 lysine 9 (H3K9me3) and, when fused to dCas9, allows us to use the programmable properties of CRISPR to target and inhibit any genomic loci of interest (Fig. 3a) [48][49][50][51] . For this experiment, each well received a different set of customdesigned guide RNAs (sgRNAs) to specifically inhibit one SE per well (i.e. arrayed CRISPRi screen) ( Fig. 3b and Supplementary Data 1). A total of 86 super-enhancers were tested plus 10 control wells. Two different sgRNAs, targeting the two highest BRD4 peak summits within each super-enhancer, were designed for each SE (see Methods) 52 . For negative controls, we used a nontargeting scrambled sgRNA in addition to a sgRNA designed to target a dormant region of the genome (Supplementary Data 1). After 72 h of epigenetic silencing, RNA was purified from each well and barcoded to specifically track which super-enhancer was probed per well (96 total barcodes). The RNA was prepped and sequenced on an Illumina platform to measure changes in gene expression as a consequence of super-enhancer inhibition (Fig. 3b).
Given our intent to survey as many super-enhancers as possible and to increase the chances of finding those that exhibited the most profound effects on gene expression, we decided to probe each SE once within the 96-well setup, prioritizing breadth over the inclusion of replicates (Fig. 3b). Therefore, a traditional differential gene expression analysis pipeline (requiring the use of replicates) had to be eschewed in favor of something better able to handle our experimental setup. We took inspiration from previous analyses performed on largescale perturbation databases, such as the Connectivity Map project (CMap) 53 , and chose to focus on relative changes in rank for each gene (uprank or downrank) rather than traditional differential gene expression analysis or absolute expression counts. The resulting changes in rank could then be investigated across the entire dataset by iterating through a series of rank change cutoffs, identifying super-enhancers that affected   significantly more genes at a particular cutoff as compared to the negative control wells (based on an empirical false discovery rate of 0.1) (see the "Methods" section). Any genes detected at these thresholds could then be tentatively assigned as target genes to each SE ( Fig. 3 and Supplementary Data 6). To investigate whether a traditional relative expression approach would have identified similar target genes, we determined the log2-fold change of every gene for each SE relative to the controls. We then assessed the relationship between gene expression determined by the relative change in rank and relative expression for each SE. In every case, the correlation between log2-fold change (LFC) and rank change (RC) was highest when comparing each SE to itself, as opposed to all other SEs on the screen, suggesting that differential gene expression calculated in both ways gave similar results ( Supplementary Fig. 3a). Notably, some of the correlations were much stronger than others, leading us to focus on SEs with an LFC versus RC correlation value above the mean. Of particular interest was super-enhancer 14 (SE14) which exhibited an LFC vs. RC correlation value of 0.95 (the highest in the entire dataset), suggesting particularly robust results for this SE ( Supplementary  Fig. 3b). Confident that our rank change approach was adequately supported by this comparison, we proceeded to look for other SEs that exhibited profound effects on gene expression. First, we focused on a number of summary analyses from the CRISPRi screen. The median number of genes downregulated by each SE was four and there were a few salient SEs that affected a much larger number of genes (Fig. 3c). Interestingly, there was only a weak correlation between the number of differentially regulated genes and SE size (Fig. 3e) or enrichment of H3K27ac (Fig. 3f), suggesting that the effects on gene expression are not merely a function of size. Of note, super-enhancer 60 (SE60) was in the bottom half in terms of size, but it affected the greatest number of genes. Therefore, we felt it prudent to understand the specificity of our CRISPRi targeting process and empirically determine the extent of spreading of the repressive H3K9me3 mark upon dCas9-KRAB binding. To this end, we performed H3K9me3 ChIP-seq in ovarian cancer cells transfected with SE60 targeting sgRNAs versus non-targeting sgRNAs. Differential binding analysis revealed that only our region of interest (SE60) was significantly enriched for the H3K9me3 signal upon transfection of the targeting sgRNAs but not the scramble nontargeting sgRNA ( Supplementary Fig. 4c, d). Additionally, there was an increase in the H3K9me3 signal at each of our two SE60 sgRNA locations, suggesting both guides successfully delivered dCas9-KRAB to the SE target sites ( Supplementary  Fig. 4a, b). We found that the region with increased H3K9me3 was about 20 kb, spreading~10 kb from each sgRNA target site, enough to cover the entire SE. There also did not appear to be an increase in signal at the computationally predicted off-target sites, suggesting that the guides for SE60 were highly specific ( Supplementary Fig. 4e, f). Taken together, these results validate our method for designing sgRNAs and reinforce that the observations from the screen (specifically for SE60) were not due to off-target effects ( Supplementary Fig. 4).
Having supported the validity of our CRISPRi assay, we next wanted to examine the patterns of gene expression that resulted from the screen. To accomplish this, we utilized two clustering methods, K-means clustering for the target genes and unsupervised-hierarchical clustering for the SEs. We found that the differentially regulated genes could be divided into three optimal clusters that represent distinct gene expression pathways in cancer cells (Fig. 3d, g). Conversely, the SEs can be divided into 10 distinct clusters with shared patterns of gene expression (Supplementary Data 1). More specifically, CRISPRi targeting of the SEs in clusters 2-4 (containing SE14 and SE60) caused decreases in the expression of genes enriched for pathways such as KRAS signaling, estrogen response (both early and late), and epithelial to mesenchymal transition (EMT). In contrast, SEs in clusters 5-10 maintain some similarities (KRAS and early estrogen response) but also have a unique role in the regulation of the JAK-STAT pathway and immune-related pathways (Fig. 3g). Taken together, the CRISPRi screen in conjunction with our CNV analyses have allowed us to comprehensively determine which SEs have the most profound effects on gene expression and inform us of the enhancers that likely regulate key gene pathways in ovarian cancer. Based on these results, two salient SEs, SE60 and SE14, were selected for follow-up experiments.
Deletion of SE60 and SE14 causes dysregulation of oncogenic gene expression pathways leading to reduced proliferation and migration of cancer cells. Perturbation of SE60 affected the greatest number of genes in the CRISPRi screen. In addition, amplification of this SE in ovarian cancer patients is prognostic of worse patient outcomes, nominating it as a super-enhancer that is associated with oncogenic processes (Fig. 4a, b). Therefore, we wanted to experimentally determine whether SE60 drives critical gene expression programs in ovarian cancer. To that end, we designed sgRNAs flanking the BRD4 peak summit of the largest constituent enhancer within SE60 and generated three independent CRISPR-Knockout (KO) clones resulting from~1700 to 1800 bp deletions ( Fig. 4c and Supplementary Fig. 5).
Global changes in gene expression resulting from each SE60 KO clone were measured using RNA-seq. Differential expression analysis using DESeq2 54 revealed pervasive changes in gene expression with 660 genes being detected as significantly downregulated and 1090 genes being upregulated at a strict confidence threshold (adjusted p-value of 0.0005) (Fig. 4d, e, and Supplementary Data 7). Pathway analysis of the top 100 significantly downregulated genes (determined by p-value) identified significant enrichment in cell cycle progression, quiescence, metastasis, differentiation, and KRAS-signaling ( Fig. 4f) [55][56][57] , further suggesting SE60 drives critical gene expression programs in ovarian cancer. This observation was supported by clinical analysis of these predicted target genes, where increased expression of the top 100 SE60 target genes is associated with worse clinical outcomes in HGSOC patients (utilizing data from 15 ovarian cancer datasets), even after considering additional clinical variables such as stage and grade (Fig. 4i, Supplementary   Fig. 3 Systematic epigenetic silencing of ovarian cancer enriched super-enhancers using CRISPRi (dCas9-KRAB) coupled to multiplexed RNA-seq. a Cartoon showing sgRNA-guided dCas9-KRAB epigenetic silencing of a SE via enrichment of the repressive histone modification H3K9me3. b Experimental setup for the CRISPRi screen in a 96-well plate (left). Western blot showing OVCAR3 cells engineered to stably express dCas9-KRAB (right). Results are representative of 3 independent experiments. dCas9-KRAB expressing OVCAR3 cells were plated in each well. One SE was targeted per well (86 SEs plus 10 control wells). After 72 h of enhancer silencing, changes in gene expression were measured using barcoded RNA-seq. Two sgRNAs were custom designed for each SE and transfected into each corresponding well. c Horizontal bar chart showing the number of downregulated genes for each SE. SE60 and SE14 were selected for further analysis as described in the text and are indicated by arrows. d K means clustering elbow plot used to determine the optimal number of gene clusters across significant DEGs for all SEs pulled from the screen analysis. The "elbow" determines the ideal cluster number which was chosen as 3. e Scatterplot comparing SE size versus the number of downregulated genes. There is no correlation between SE size and the number of target genes. SE60 and SE14 were selected for further analysis as described in the text and are indicated by arrows. f Scatterplot comparing H3K27ac enrichment versus the number of downregulated genes. There is not a strong correlation between H3K27ac enrichment and the number of target genes. g Heatmap representing the unsupervised hierarchical clustering of all SEs (clusters 1-10 under the dendrogram) and controls in the screen across all screen DEGs (left). The boxes on the right denote the three K-means clusters. MSigDB pathway analysis describes the functions the genes in these clusters are involved in (right). SE60, SE14, and the two negative controls are denoted at the bottom of the plot. Source data are provided as a Source Data file.
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-31919-8 ARTICLE Fig. 6, and Supplementary Data 10). The notion that SE60 plays a key role in ovarian cancer was further validated by the effects that deletion of this SE had on cancer cell proliferation and migration (Fig. 4g, h).
To substantiate our approach for identifying clinically relevant super-enhancers, we selected an additional candidate from the CRISPRi screen for validation. SE14 was chosen because (1) it has the highest correlation between LFC and RC differential gene expression analysis from the CRISPRi screen, (2) it is in the top four SEs that affected the greatest number of genes, and (3) its amplification portends a worse clinical outcome in ovarian cancer patients (Fig. 5a, b). To investigate the functional role of SE14, we . Taken together, these results validate our approach to identifying clinically relevant SEs and highlight the importance of these two SEs in ovarian cancer. Next, we investigated whether the mechanistic roles of SE60 and SE14 on cell proliferation and migration were due to direct or indirect target gene regulation.
3D-chromatin interactions defined by Hi-C in ovarian cancer cells establish direct target genes for SE60 and SE14. The significant effects on proliferation and migration caused by CRISPRbased deletion of SE60 and SE14 led us to investigate if these biological phenotypes were caused by direct or indirect regulation of target genes. We reasoned that direct target genes would exhibit increased chromatin looping interactions with the SE, whereas indirect target genes would be downstream of an effector gene that was directly regulated by the SE. To enable unbiased measurement of interaction frequencies between each superenhancer and its target genes, we performed Hi-C in OVCAR3 cells to comprehensively annotate chromatin interactions across the ovarian cancer genome 58 . In order to maximize the breadth of this analysis, we focused on the target gene set detected from the CRISPR-KO experiments that represented the most statistically robust gene set for each SE, resulting from 4 replicates of RNAseq across three independent knockout clones for each SE. Moreover, the constitutive perturbation of each SE, caused by CRISPR-based deletion, gave rise to consistent gene expression patterns that resulted in marked biological phenotypes, thus facilitating integration with the Hi-C data.
Since Hi-C is highly dependent on distance, we limited our search space to genes located on the same chromosome as the super-enhancers (cis genes) in order to get an accurate metric of interaction frequency 59,60 . To perform this analysis, we quantified the interaction frequency between each SE and its downregulated genes upon SE deletion. This was compared to a control dataset consisting of 100 permutations of distance-matched gene sets that exhibited no significant changes in gene expression upon SE deletion (see the "Methods section). This enabled us to compare distributions of interaction frequency measurements between each SE and a random set of genes based entirely on genomic distance. Direct targets were defined as SE-gene pairs with an observed/ expected contact frequency greater than the 75th percentile of the control/background distribution. Overall, we observed that the target genes for each SE had a higher interaction frequency with their cognate SE compared to distance-matched genes found on the same chromosome (Fig. 6a, f).
We identified one cis direct target gene and four cis indirect target genes for SE60 (Fig. 6b). Of note, the cis direct target gene for SE60, RAE1, has previously been associated with invasion in ovarian cancer and has been shown to promote EMT in breast cancer 61 . In addition, increased expression of this gene portends a worse outcome in HGSOC patients 34,35 Fig. 4 CRISPR-knockout of super-enhancer 60 leads to profound changes in gene expression and reduced proliferation of cancer cells. a Genome browser view of SE60 (dashed red box) and the surrounding region showing enrichment of BRD4, H3K27ac, and ENCODE H3K27ac signal. b Kaplan-Meier plots of copy number amplification over each SE60 overlapping 15 kb windows versus disease-specific survival in TCGA HGSOC patients. Significance was assessed using a log-rank test and Cox proportional hazards model. c Top: Cartoon showing CRISPR-knockout (KO) of SE60. Bottom: Genotyping PCR showing successful heterozygous knockouts of SE60-representative of three independent experiments. d Unsupervised hierarchical clustering heatmap of all 1750 significant DEGs (Pearson correlation, adjusted p-value > 0.0005 at any fold change) between wild-type and SE60 KO cells measured by RNA-seq. e PCA plot showing the variance landscape of WT and KO samples. f Pathway analysis using CancerSEA and MSigDB of the 100 most significant DEGs detected (as determined by Kolmogorov-Smirnov style GSVA followed by Spearman's rank test with BH FDR correction). The red line denotes the metric for a p-value of 0.05 converted into the −log10 scale. g Proliferation assays of three independent SE60 KO clones (represented in the RNA-seq data) versus wild-type OVCAR3 cells. Results are shown as fold change compared to day 2 values. The statistically significant differences (determined by a twosided Student's t-test) are provided in red text. n = 4 biological replicates. Data are shown as mean ± SEM. h Cell Migration assays of three independent SE60 KO clones versus wild-type OVCAR3 cells. Microscope brightfield images of the growth after 24 h (left). Bar chart representation of cell count after 24 h, statistically significant differences (as determined by a two-sided Student's t-test) are provided in red text (right). Data are shown as mean ± SEM. n = 4 biological replicates. i Kaplan-Meier plot 34,35 showing the clinical significance of the top 100 downregulated genes after SE60 KO. Significance was assessed using a log-rank test, significant p-values are denoted in red text. Source data are provided as a Source Data file.  (Fig. 6c, d, Supplementary Fig. 6, and Supplementary Data 10). Notably, RAE1 was also predicted as a target of SE60 by the CNVeQTL analysis ( Supplementary Fig. 9a, b, and Supplementary Data 8). When looking at Hi-C contact frequency across chromosome 20, we notice a marked increase in contact between the RAE1 locus and the SE60 locus as compared to the background (Fig. 6e). This suggests that the decrease in migration detected upon SE60 deletion is due, in part, to its direct regulation of RAE1. We suspect that there may exist more direct target genes for SE60 located on other chromosomes based on the reference genome, that may be translocated to chromosome 20 in the ovarian cancer genome. However, these trans-chromosomal interaction frequencies are technically more challenging to detect via Hi-C.
Interestingly, we identified a much greater number of cis direct target genes (28 genes) and cis indirect targets (62 genes) for SE14 ( Fig. 6g and Supplementary Data 8). Pathway analysis of the cis direct targets revealed key roles in cell cycle progression,  quiescence, invasion, differentiation, metastasis, and stemness (Fig. 6h). Kaplan-Meier analysis of this gene signature highlighted a statistically significant decrease in survival for patients that had high expression of these genes (Fig. 6i, Supplementary  Fig. 6, and Supplementary Data 10). Likewise, 8 of these cis direct targets had been predicted from our CNVeQTL analysis reinforcing the utility of CNVeQTLs to predict cis-direct targets ( Supplementary Fig. 9c, d). Through our analysis of these of cis direct targets, we identified examples of both close-range (EPHA2) and distant (MAB21L3) connections to SE14 ( Fig. 6j and Supplementary Fig. 10). Interestingly, there were genes within very close proximity to SE14 (such as ARHGEF19) that showed no evidence of interaction or differential gene expression. These data implicate SE14 as being directly involved in both proliferation and migration, as well as other key processes in ovarian cancer. To validate our Hi-C analysis, we performed an orthogonal method, called Activity-by-Contact model (ABC) 62 , to identify interactions across the genome by integrating Hi-C data with measures of chromatin activity. This analysis allowed us to determine the top interacting cis-genes for both SEs regardless of the distance between both features. The results of the ABC analysis recapitulate those from our initial Hi-C analysis and independently nominated RAE1, EPHA2, and MAB21L3 as the genes most likely to interact with SE60 and SE14, respectively (Supplementary Fig. 11 and Supplementary Data 9). This analysis highlighted our ability to delineate true targets and suggests that both SE60 and SE14 are directly involved in clinically relevant processes in ovarian cancer.
SE60 and SE14 are specifically active within the epithelial cancer cell fraction of human HGSOC tumors as revealed by single cell genomics. Our previous experiments had demonstrated that these SEs are preferentially amplified in ovarian cancer patients and that they regulate gene expression pathways that govern the proliferation and migration of cancer cells. As a final validation experiment, we wanted to determine if SE60 and SE14 were specifically active within the cancer cell compartment of human HGSOC tumors and if their target genes are also active within the same cell type. To test this, we analyzed matched single-cell RNA-seq (scRNA-seq) and single-cell ATAC-seq (scATAC-seq) data from two HGSOC patients previously generated in our lab ( Supplementary Fig. 12) 43 . We annotated seven distinct cell types present in these tumors by both scRNA-seq and scATAC-seq and identified the cancer cell population using the FDA-approved biomarker CA125 (also known as MUC16) (Fig. 7a, b) 63 . We found significant enrichment of RAE1, a SE60 cis direct target, and EPHA2, a SE14 cis direct target, within the cancer cell fraction as compared to the normal cell fraction (Wilcoxon Rank Sum tests, Bonferroni-corrected p-values < 2.2e −308 and average logFC ≥ 0.1) (Fig. 7b).
In order to assess whether SE60 and SE14 are uniquely active in ovarian cancer, we next leveraged the scATAC-seq data. These data showed significantly increased chromatin accessibility at three constituent enhancers of both SE60 and SE14, specifically within the cancer epithelial cell fraction as compared to the stromal compartments of these tumors (Wilcoxon Rank Sum tests, Benjamini-Hochberg FDR ≤ 0.10 and Log2FC ≥ 0.25) (Fig. 7c). Additionally, both HGSOC patients showed this pattern, suggesting that activation of these SEs is a common feature of HGSOC biology. While there is previous evidence from ENCODE that these regions contain regulatory elements in normal epithelial tissue, it appears that there is significantly more accessibility of these super-enhancers in ovarian cancer cells (Fig. 7c).
In order to investigate what transcription factors might be involved with these super-enhancers, we performed motif enrichment analysis using FIMO sequence analysis 64 . To provide confidence to the TF motif calls, we investigated the expression of the predicted TFs within the cancer epithelial cells. Within the cancer-enriched constituent enhancers of SE60, we found SOX4, ATF4, and YY1 as the top three predicted transcription factors. Of note, YY1 is known as an integral component of enhancer-promoter loop interactions and is a hallmark of an active enhancer 65 . Similarly, we detected binding motifs for ELF3, KLF, and JUN in the cancer-enriched constituent enhancers of SE14. Notably, ELF3 has been previously associated with vascular inflammation, tumorigeneses, epithelial differentiation, and the ERRB3 pathway providing additional evidence of the importance of this SE (Fig. 7d, e) 66,67 . Taken together, these data suggest that these SEs and their target genes show more activity in the cancer cell fraction of HGSOC tumors and serve to validate our computational pipeline for the identification of clinically relevant super-enhancers.

Discussion
Every year, an estimated 22,000 new cases of ovarian cancer will be diagnosed and around 14,000 women will die as a result of this disease 1 . The paucity of known drivers for ovarian cancer makes identifying at-risk individuals very difficult and has led to a lack of effective targeted therapies. Thus, platinum-based chemotherapy coupled with surgery remains the standard of care 68 . Given their critical functions in controlling gene regulation, enhancers are often required to achieve the levels of transcriptional activity needed to sustain cancer cells and have been shown to play an integral part in cancer development and patient survival. Additionally, super-enhancers have demonstrated the capacity to regulate many critical pathways for the development and maintenance of the cancer cell state as well as influence therapeutic resistance [17][18][19][20] .
With the advent of therapeutics designed to inhibit various epigenetic factors that convey functionality to enhancers, it is now possible to exploit the dependency of cancer cells on transcription as an effective strategy for treating therapeutically recalcitrant cancers such as ovarian cancer 27,69 . For example, the Bromodomain and Extra-Terminal motif inhibitors (BET inhibitors; such as JQ1) designed to interfere with the functions of bromodomaincontaining proteins like BRD4 have shown promise in several pre-clinical models of cancer, although their efficacy in a clinical setting is still unknown 25,70 . Nonetheless, investigating enhancers with high BRD4 enrichment can lead to the identification of biomarkers, druggable targets, and an improved understanding of ovarian cancer. Notably, the expression of BRD4 is highest in ovarian cancer as compared to every other cancer type represented in The Cancer Genome Atlas and high expression portends a worse outcome in ovarian cancer patients (Fig. 1). Thus, we reasoned that co-enrichment of BRD4 and H3K27ac can be used as a surrogate to find SEs driving oncogenic processes in ovarian cancer. This was substantiated by the observation that the SEs identified in our study were preferentially copy number amplified in ovarian cancer patients and that some amplification events were themselves predictive of worse clinical outcomes (Fig. 2). Additionally, our CNVeQTL analyses across HGSOC patients demonstrate that the activity of these super-enhancers is pervasive. This is perhaps not surprising since genomic instability is a hallmark of ovarian cancer and several studies have demonstrated that somatic mutations at specific regulatory elements in the ovarian cancer genome play a pivotal role in subtype determination and overall progression 2,71 . Furthermore, the dysregulation of genomic architecture in ovarian cancer may allow for cancer cells to hijack existing enhancers for oncogenic purposes. In fact, several examples of enhancer hijacking exist in other types of cancer such as Burkett's Lymphoma, B-Cell Lymphoma, and Glioblastoma 9,72,73 . Overall, these findings suggested that a number of our identified SEs were amplified for biologically meaningful reasons. Rather than limiting our study to the standard taxonomic listing of super-enhancers, we used three orthogonal approaches to define the regulatory logic of SEs in ovarian cancer-CRISPRi, CRISPR-KO, and Hi-C. The CRISPRi screen allowed us to systematically determine the target genes for each of the top 86 most active SEs (Fig. 3). While most CRISPR screens involve a pool of sgRNAs and rely on a cellular endpoint (such as proliferation) to be able to capture the relative abundances of remaining sgRNAs, our screen was customized to provide a readout of gene expression for each super-enhancer. We knew, a priori, which sgRNAs were used and which SEs were affected in each well. On average, we found that each SE perturbation resulted in the downregulation of about four genes and the total number of genes regulated by each SE was not a function of size or enrichment of H3K27ac or BRD4. In fact, SE60 was in the bottom quartile of super-enhancers in terms of size and H3K27ac enrichment, but it had the most profound effects on gene expression. Therefore, we reasoned that SE60 harbored the most potential for further study due to its likely role in regulating genes that contribute to the pathology of ovarian cancer. While the goal of the CRISPRi screen was to broadly investigate the effects on gene expression across a large cohort of super-enhancers, we recognize that the CRISPRi screen was underpowered to definitively establish target genes for each SE. Thus, we elected to perform CRISPR-KOs of SE60 and SE14 to enable robust target gene detection.
CRISPR-KO of SE60 and SE14 had dramatic effects on gene expression programs involved in Quiescence, Metastasis, and Invasion, among other important pathways. Moreover, the gene sets for both SEs were associated with poor outcomes in HGSOC patients. This was supported by both proliferation and migration defects in the SE60 and SE14 knockout cells (Figs. 4 and 5). We note that there were hundreds of genes differentially regulated upon deletion of these two SEs, and that there was a modest overlap with the differentially expressed genes detected via CRISPRi. This is perhaps due to the technical nuances of CRISPR-KO (a constitutive genomic deletion arising from an individual clone) versus CRISPRi (a transient epigenetic inhibition of the enhancer locus) that may affect the target genes identified for a particular SE. In fact, the field as a whole has wrestled with the best way to assign target genes to enhancers, especially considering the genomic rearrangements observed in cancer cells. Thus, we reasoned that direct chromatin interactions between the SEs and their target genes (as measured by Hi-C) would give confidence to the annotation of target genes (Fig. 6).
Overall, the downregulated genes upon SE deletion showed higher interaction, both nearby and across long distances, with the SE as compared to the distance-matched control gene set. Importantly, several cis direct target genes are involved in oncogenic pathways and perhaps could serve as prognostic indicators or biomarkers in the future (Fig. 6, Supplementary  Fig. 6, and Supplementary Data 10). We note that there may exist more target genes for each SE located on other chromosomes, however, we found no distinguishable interactions (via Hi-C) between our SEs and their target genes located on different chromosomes. This suggests that genes found on different chromosomes are likely indirect target genes where the SE directly regulates a gene in cis that is co-linear with the trans target gene. In the absence of Hi-C data for determining direct target genes, we posit that evidence from two orthogonal experiments (such as CRISPRi and CRISPR-KO or inclusion of reporter-based enhancer assays) would yield high confidence results since genes detected by multiple assays are agnostic to the technical nuances of each. In fact, a logical framework to describe the level of support needed to definitively annotate an enhancer and its bona fide target genes have been recently proposed 74 , and its implementation would yield a catalog of enhancers with confidently linked target genes.  6 Hi-C analysis detects direct targets of SE60 and SE14 supporting direct roles in invasion, differentiation, and metastasis. a Distribution of Hi-C counts (contact frequency) between SE60 downregulated cis-genes and the SE60 locus (left) (n = 4), blue points/genes are direct targets. Distribution of Hi-C counts of a background set of 100 distance-matched gene sets (right) (n = 400). The dashed line denotes the cutoff for direct target genes in the experimental sample. b Table displaying the number of direct and indirect cis-target genes of SE60 determined by Hi-C. c CancerSEA analysis of the direct target (p-values determined by GSVA followed by Spearman's FDR corrected rank test). The red line denotes a p-value of 0.05. d Kaplan-Meier plot of RAE1 34,35 . The red line denotes patients with high expression of RAE1, significant p-values are denoted in red. e Hi-C contact heatmap showing the interaction between RAE1 and the SE60 locus (red square). The colored gene names correspond to the fold change of significantly affected genes after CRISPR-KO. f Distribution of all Hi-C counts between SE14 downregulated cis-genes and the SE14 locus (left) (n = 90), blue points/genes are direct targets. Distribution of Hi-C counts of a background set of 100 distance-matched control gene sets (right) (n = 9000). The dashed line denotes the cutoff for direct target genes. g Table displaying  Finally, both SE60 and SE14 were found to have a statistically significant increase in chromatin accessibility within the cancer cell fraction of human HGSOC tumors at single-cell resolution (Fig. 7), further suggesting that the SEs that we identified are not merely cell line specific. This validates our enhancer identification pipeline and reveals that certain super-enhancers involved in growth and migration are preferentially enriched and amplified in cancer cells. In addition, we found that cis direct target genes annotated for each SE (such as RAE1 and EPHA2) were more highly expressed in the cancer cells compared to the stromal/nonmalignant cells within HGSOC tumors. Collectively, these results expound the idea that super-enhancers themselves and the genes  Engineering dCas9-KRAB expressing OVCAR3 cells. Lentivirus was packaged in HEK-293T cells and contained the Lenti-dCas9-KRAB-blast vector 75 (Addgene plasmid #89567). Cells seeded in a T75 flask were transfected with the following: 6.67 µg Lenti-dCas9-KRAB-blast, 5 µg psPAX2 (Addgene, 12260), and 3.33 µg PMD2G (Addgene, 12259) using Fugene 6 (Promega, E2691) following the manufacturer's protocol. 48-72 h post-transfection, the lentivirus containing supernatant was harvested and concentrated using the Lenti-X Concentrator (Takara, 631231) following the manufacturer's protocol. To transduce OVCAR3 cells, cells were seeded at 50,000 cells/well in a six-well plate and treated with the harvested lentivirus in RPMI media with 10% FBS and 10 µg/mL polybrene (Millipore, TR1003G). After 72 h, transduced cells were placed in RPMI selection media with 3 µg/mL blasticidin for 7 days. After batch selection, OVCAR3-KRAB cells were collected for Western Blot to validate the presence of dCas9-KRAB. Cells were lysed with following lysis buffer: 50 mM Tris-HCl (pH 8), 0.5 M NaCl, 1% NP-40, 0.5% sodium deoxycholate, 0.1% SDS and 1× protease inhibitor. The β-tubulin antibody (Abcam, ab6046) was diluted at 1:5000 in 5% BSA in TBST and incubated overnight at 4°C. The Cas9 antibody (7A9-3A3) (Santa Cruz, sc-517386) was diluted at 1:1500 in 5% BSA in TBST and incubated overnight at 4°C. The secondary antibodies (Donkey anti-rabbit IgG HRP-linked (GE, NA934) and Donkey anti-mouse IgG HRP-linked (Invitrogen, PA1-28748)) were diluted at 1:5000 in 5% BSA in TBST.

TF name Motif ID Strand p-value q-value
CRISPRi screen sgRNA design. For sgRNAs targeting super-enhancers, target regions were chosen by selecting two regions within the super-enhancer with the highest BRD4 enrichment and clear H3K27ac signal. For each super-enhancer region, sgRNAs were designed using the CRISPOR web tool 52 taking into account the specificity and off-target scores. If all suggested sgRNA sequences to a region had low specificity scores, a second sgRNA was instead designed to target the third highest BRD4 peak. Two sgRNAs were designed per super-enhancer to be transfected together. Genomic coordinates for all super-enhancers and their sgRNA sequences are found in Supplementary Data 1. sgRNA oligos were ordered from Integrated DNA Technologies. The negative control sgRNAs (Scramble1 and Scramble2) were previously published 76 .
Additionally, the sgRNA stem-loop was extended and an A-U base pair flip was utilized to improve sgRNA stability and assembly with dCas9 78 . The vector cloning protocol was adapted from Feng Zheng's group 79  Super-enhancer knockout with CRISPR-Cas9. Deletion of super-enhancers was performed using CRISPR-cas9 following existing protocols [79][80][81] . In short, the BRD4 peak summit was targeted for each super-enhancer and sgRNA target sites were selected using the CRISPOR web tool 52 . The guide RNA sequences and their genomic coordinates can be found in Supplementary Table 1 Table 2). Deletion of SE14 resulted in a~2500-2800 bp deletion. Deletion of SE60 resulted in a~1700-1800 bp deletion.  For SE60KO3, SE14KO2, SE14KO3, scramble1-KRAB, and SE60-KRAB, libraries were created and sequenced by Novogene. These libraries underwent 150 bp paired-end sequencing on an Illumina NovaSeq 6000 instrument.
ChIP-seq. OVCAR3-KRAB cells were transfected with sgRNAs targeting either scramble1 (non-targeting) or SE60 (2 pooled sgRNAs) following the same protocol mentioned above for "CRISPRi for SE14 and SE60." For each of the two replicates conducted per condition, 1-2 million cells were used for fixation with 11% formaldehyde following Active Motif's Epigenetic Services ChIP Fixation Protocol. ChIP-seq for H3K9me3 was performed by Active Motif using H3K9me3 antibody (Active Motif, 39161) with spike-in Drosophila normalization. ChIP-seq libraries underwent 75 bp single-end sequencing on an Illumina NextSeq 5000 instrument by Active Motif.
Cell proliferation assay. Cell collections were performed at Days 2, 4, and 6. On the day of collection, cells were fixed with 10% formaldehyde and stained using a 0.1% crystal violet solution. Incorporated crystal violet was extracted using 10% glacial acetic acid and the absorbance was read at 595 nm. This procedure was conducted four times to ensure reproducibility. Results are shown as the mean fold change of Day 4 and Day 6 OD 595 nm readings compared to Day 2 OD 595 nm readings ± SEM. Statistical analysis was conducted in R using a two-sided Student's t-test.
Cell migration assay. OVCAR3 WT and SEKO cells in serum-free RPMI media were seeded to the upper chamber of a transwell insert at 60,000 cells per insert. The lower chamber contained RPMI with 10% FBS. Cells were incubated for 24 h, then all non-migrated cells were removed from the upper membrane. Cells were fixed and stained using the Hema 3 Staining Kit (Fisher Scientific, 122-911). Ten brightfield images were taken per insert and images were analyzed using the CellProfiler 4.2.1 software to count the number of cells per transwell insert. This procedure was conducted four times to ensure reproducibility. Results are shown as the mean cell count per insert ±SEM. Statistical analysis was conducted in R using a t-test.
General program versions. Unless specified, these are the versions used for scripting/analysis in R (R: 4.0.0) and Python (Python: 3.6.5) throughout the project for the bulk data analysis of CRISPRi, CRISPR-KO, CNV, and H3K27ac/BRD4 ChIP-Seq data. Unless otherwise stated, all "overlap" analysis visualization was performed using Intervene Intervene: (0.6.5) 82 .
RNA Seq: CRISPRi screen General metrics. RNA-seq was performed following the pipeline put forth by LEXOGEN in the 3' mRNA-Seq package; namely using STAR, HTSEQ, and DESEQ2. These processes will be explained in more detail below.
QC. Quality control was performed using the FastQC (version v0.11.7) tool and the results were analyzed (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). All of the metrics returned as acceptable with no clear failures. We thus proceeded with processing and analysis.
File formatting. The bam files from STAR were then indexed and sorted using functions in the SAMTOOLS package (version 1.9), namely samtools sort and samtools index 84 .
Quantification. The sorted and indexed bam files were quantified using htseq (version 0.11.2) and the gencode v29 primary assembly as a reference and with the following parameters -m intersection-nonempty, -s yes, -f bam, -r pos 85 .
Read distributions. The package RSeQC (version 3.0.0) was used to assess the distribution of reads across the genome. Specifically, the python program read_distribution.py was used with default parameterizations to create a summary of this information 86 .
Review QC. All of the alignment, counting, and cleaning program outputs were assessed with MultiQC (version 1.9) for potential issues, of which none were determined 87 . Default parameters were used.
Normalization. The count data were first normalized by removing all of the low count genes (genes with <1 count in every sample); this data was then read into DESEQ2 (DESeq2_1.30.1) 54 . Within DESEQ2 normalized by scaling and size factors followed by a VST transformation. Batch effects were addressed by utilizing the SVT program (part of the DESEQ2 package using sva_3.38.0) and variation from two surrogate variables was removed for the final analysis. The process used for this step of the analysis can be followed within the script Screen_Preprocessing.R.
Determination of DEGs. Differential gene expression was determined by utilizing a rank-based approach similar to the ranking method used by CMAP for their single replicate screens 53 . Genes were ranked in order of expression (rank 1 being the highest expressed, n being the lowest) within every sample, then all samples were aggregated and a global rank was assigned for every gene. Next, the change in rank was determined between the within-sample rank and the global rank for every gene in every sample. These changes in rank were used to build a distribution of all rank changes for eFDR analysis. The process used for this step of the analysis can be followed within the script OVCAR3_Screen_Analysis_with_Plotting_LFC_Comparison.ipynb.
Empirical false discovery rate analysis. Empirical false discovery rate, an empirically derived variation of the false discovery rate, was determined by choosing a rank change threshold and assessing the median number of genes across controls beyond that threshold as compared to a given sample 47 . For example, if there is a median of 4 genes in the controls and 40 genes in Sample A; the eFDR for this comparison would be 4/40 or 10%. The process used for this step of the analysis can be followed within the script OVCAR3_Screen_Analysis_with_Plotting_LFC_Comparison.ipynb.
Relative expression correlation analysis. A log2-fold change was calculated between all genes in a sample and the median of the controls. All genes determined as significant by the rank-based analysis (across all super-enhancers) were aggregated into one pool of genes. This pool of genes was then used to compare RC to LFC values within each super-enhancer to determine the correlation of these sets of values. The process used for this step of the analysis can be followed within the script OVCAR3_Screen_Analysis_with_Plotting_LFC_Comparison.ipynb.
Clustering. KMeans clustering analysis was used to cluster the differentially ranked gene list. Three clusters were determined as optimal by analysis of the elbow plot and these clusters were then applied to the data. Unsupervised hierarchical clustering was then used to determine the super-enhancer relationships. The process used for this step of the analysis can be followed within the script OVCAR3_ Screen_Analysis_with_Plotting_LFC_Comparison.ipynb.
Pathway analysis. Genes detected from the differential expression analysis were analyzed using CancerSEA and the molecular signatures database [55][56][57] . This program performs pathway analysis using cell-type specific information relevant to cancer based on available single-cell datasets. All of the genes in a given KMeans cluster were fed into this set of programs as a gene list and results were retrieved.
RNA Seq: CRISPR KO General metrics. RNA-Seq was performed following a similar pipeline to that used in the screen analysis with parameters adjusted to account for differences in the data (paired-end with greater depth); namely using STAR, HTSEQ, and DESEQ2. This will be expounded in more detail below.
Trimming. No trimming was needed or performed.
File formatting. The bam files from STAR were then indexed and sorted using functions in the SAMTOOLS package (version 1.9), namely samtools sort and samtools index 84 .
Read distributions. The package RSeQC (version 3.0.0) was used to assess the distribution of reads across the genome. Specifically, the python program read_distribution.py was used with default parameterizations to create a summary of this information 86 .
Review QC. All of the alignment, counting, and cleaning program outputs were assessed with MultiQC (version 1.9) for potential issues 87 . Default parameters were used and all of the reports were good.
Normalization (batch effect detection). The count data were first normalized by removing all of the low count genes (genes with <1 count in every sample); this data was then read into DESEQ2 (DESeq2_1.30.1) 54 . Within the DESEQ framework, the counts data were adjusted for scaling and size factors followed by a VST transformation. Batch effects were addressed by utilizing the SVT program (sva_3.38.0) and variation from one surrogate variable was accounted for in the DESEQ2 model. This process was completed and can be followed using DESEQ2_2021Reps_RNA_SVA_Plotting_V2.Rmd.
Normalization. The pre-VST data was used for standard in-program normalization by DESEQ2 during the differential expression analysis procedure. This process was completed and can be followed using DESEQ2_2021Reps_RNA_SVA_Plotting_ V2.Rmd.
Determination of DEGs. Differential gene expression was determined by utilizing DESEQ2 and default parameters. Genes called as differentially expressed at an FDR-adjusted p-value less than 0.0005 were identified and collected for analysis and figure making. This process was completed and can be followed using DESEQ2_2021Reps_RNA_SVA_Plotting_V2.Rmd.
Pathway analysis. Genes detected from the differential expression analysis were analyzed using CancerSEA and the molecular signatures database [55][56][57] . This program performs pathway analysis using cell-type-specific information relevant to cancer based on available single-cell datasets. The top 100 most significant downregulated genes from differential expression analysis were fed into this program as a gene list and results relevant to ovarian cancer were retrieved. This process was completed and can be followed using DESEQ2_2021Reps_RNA_SVA_Plotting_V2.Rmd.
Survival analysis. To perform survival analyses we made use of the KM plotter tool 34,35 . This tool allows a user to look at the effect that the expression of induvial genes or a gene set has on overall survival across a number of cancer patients from 15 ovarian cancer datasets. We looked at the top 100 genes ordered by adjusted Pvalue (the top 100 most significant genes) as a set (using the median expression of the whole group); and/or looked at genes individually.
CRISPRi RNA-Seq analysis General metrics. RNA-Seq was performed following a similar pipeline to that used in the screen analysis with parameters adjusted to account for differences in the data (paired-end with greater depth); namely using STAR, HTSEQ, and DESEQ2. This will be expounded in more detail below.
Trimming. No trimming was needed or performed.
File formatting. The bam files from STAR were then indexed and sorted using functions in the SAMTOOLS (version 1.9) package, namely samtools sort and samtools index 84 .
Read distributions. The package RSeQC (version 3.0.0) was used to assess the distribution of reads across the genome. Specifically, the python program read_distribution.py was used with default parameterizations to create a summary of this information 86 .
Review QC. All of the alignment, counting, and cleaning program outputs were assessed with MultiQC (version 1.9) for potential issues 87 . Default parameters were used and all of the reports were good.
Normalization. The pre-VST data was used for standard in-program normalization by DESEQ2 during the differential expression analysis procedure. This process can be followed using DESEQ2_RNA_Plotting_CRISPRi_Analysis_Revised.Rmd.
Determination of DEGs. Differential gene expression was determined by utilizing DESEQ2 and default parameters. Genes called as differentially expressed at an FDR adjusted p-value less than 0.0005 were identified and collected for analysis and figure making. This process can be followed using DESEQ2_RNA_Plotting_CRISPRi_ Analysis_Revised.Rmd.
Pathway analysis. Genes detected from the differential expression analysis were analyzed using CancerSEA and the molecular signatures database [55][56][57] .
Survival analysis. To perform survival analyses we made use of the KM plotter tool 34,35 . This tool allows a user to look at the effect that the expression of induvial genes or a gene set has on overall survival across a number of cancer patients from 15 ovarian cancer datasets. We looked at the top 100 genes ordered by adjusted Pvalue (the top 100 most significant genes) as a set (using the median expression of the whole group); and/or looked at genes individually.

Copy number analysis
Gathering. The copy number and RNA-seq data for this analysis were downloaded from the TCGA repository Firebrowse (http://firebrowse.org/) which contains the data used in the TCGA analysis of ovarian cancer 10 . We used the TCGA patient barcodes to determine if a tumor was from normal tissue or cancer patients. Samples were subset based on these barcodes to select for tumors. Additionally, for the CNVeQTL analysis, samples unique to each dataset (RNA or Copy Number) were removed. To perform this, we looked for matching patient identifiers between RNA-seq and copy number data and kept any data with ID overlaps.
Windowing. The autosomal (Chr 1-22) genome (hg19) was divided into 15 kb bins using python. We decided to use a sliding window size of 15 kb based on the overall size distribution of our super-enhancers. Since the median size of our superenhancers is 21 kb, we wanted a window size similar to the median size but smaller, as smaller windows allow for better resolution. We settled on 15 kb as being close to the median size and small enough to give us good resolution, yet large enough to be computationally feasible (smaller window sizes create larger datasets and increase the computational burden of assigning signals and analyzing the data). This process was completed using Split_Genome_into_windows.ipynb.
Super enhancer overlap. Bedtools (version v2.25.0) intersect (one bp overlap) was used to create a subset of the whole genome 15 kb sliding windows that overlapped the super-enhancer regions. This gave us two data sets, one being the whole genome sliding windows and the other being SE overlapping sliding windows (a subset of the whole genome group).
Copy number assignment-whole genome. Patient copy number was assigned to each 15 kb window for every chromosome individually using the script OVLP_CNV_Whole_Genome.py. If the patient data overlapped a sliding window by at least one base pair, signal from the patient was assigned to this window. Once this was performed for every chromosome individually, the chromosome data were aggregated using Combine_CNV_Chr_Files.ipynb.
Copy number assignment-super-enhancer overlap. Patient copy number was assigned to each 15 kb bin for every chromosome individually using the script SEOVLP_CNV.py. If the patient data overlapped a sliding window by at least one base pair, signals from the patient was assigned to this window. Chromosome data was aggregated using Combine_CNV_Chr_Files.ipynb.
CNVeQTL analysis. Copy Number Expression QTL was identified using MatrixQTL where the SE overlapping CNV windows were defined as the "SNPs" and the matching RNA-Seq data served as the Expression dataset 46 . Of note, genes with over 100 NA, missing, or 0 values were removed from this dataset prior to analysis. The CNV and RNA data were also converted into float values for ease of use in MatrixEQTL. CNVeQTL was identified using the linear MatrixEQTL algorithm on the original data with a P-value threshold of 1e−3. In order to determine significance, the null hypothesis was induced and used to determine an empirical FDR. The null hypothesis, in which there is no association between specific copy number regions and gene expression, was induced by randomly permuting the column assignments of the RNA-seq data, the CNV data was left alone. This maintains the variance structure of the CNV data and merely changes which CNV data column gets matched with a given RNA-Data column. For example, CNV columns 1-4 (corresponding to patients 1-4) might now be matched with RNA columns 30, 75, 6, and 210; this allows us to use the same overall data and investigate what happens where there is no link between CNV and RNA values (as patient 1's CNV values should be random in relation to the gene expression of patient 30). MatrixEQTL was then run using the original copy number data and the new column shuffled RNA-seq data; this shuffling and running of MatrixEQTL were performed 100,000 times. The median number of significant eQTLs detected across all 100k null conditions was used as the numerator for the empirical false discovery rate analysis, with the experimental results being the denominator. There is some variability in eFDR, as no seed was set and the permutations are random, but all repeats of 100k (3 repeats or 300k trials) returned an eFDR <0.1 or 10%. This process was completed using CNV_eQTL.R.
Determining super enhancer amplification. In order to assess whether the superenhancer regions were amplified, we compared the distribution of CNV values in the super-enhancer overlapping sliding windows with the whole genome by subsetting and direct comparison. We performed 10k random subset comparisons, and one direct comparison. In any given comparison, we took the 336 superenhancer overlapping windows and then randomly drew 336 windows from the whole genome background; these two sets were then compared for significant differences using a Welch's one-sided t-test. This analysis allowed us to determine if the super-enhancer overlapping group was significantly amplified relative to the randomly drawn subset. For the direct comparison, we took all 336 SE overlapping windows and directly compared the CNV values across these windows to thẽ 192,000 total regions using the same t-test metric. This process was completed using the script OVCAR_CNV_Comparison_Final.R.
Survival analysis. The effect of amplification of these regions on overall survival in patients was calculated using the Kaplan-Meier log rank change test and the Cox proportional hazards model 32 . The survival data were downloaded from the TCGA and the patient ID was mapped back to the CNV values for each patient 44 . These datasets were then combined into a single set formatted as described in CNV_KM_Plots.R. This combined survival and copy number dataset was then analyzed using the functions built in CNV_KM_Plots.R to provide a metric of significance for each 15 kb copy number region.
ChIP Seq (OVCAR3 BRD4 and H3K27ac) Data acquisition. Publicly available ChIP-Seq data were downloaded from the SRA database associated with GSE101408 (experimental OVCAR3 H3K27ac condition) using fastq dump (version 2.9.0) 30 . This process was repeated to get BRD4 binding data for DMSO-treated OVCAR3 cells as well as the input control from GSE77568 28 .
Processing. The following steps were used to process each file separately (H3K27ac, BRD4 ChIP, and BRD4 sample input). At the peak calling step, the BRD4 ChIP data was informed by the processed input control. As there was no input provided for the H3K27ac data, no input was processed or utilized for this sample. Alignment. The fastq files were aligned to hg19 using Bowtie2 (version 2.3.5.1) with default parameters 89 . The output sam files were then converted to bam files using samtools and sorted/indexed.
Processing Bam files (removing duplicates). The duplicate reads marked by Picard MarkDuplicates 2.11.0 were then removed by samtools using the command samtools view -F 1804 -b in.bam > clean.bam.
Create tagAlign files. A tagAlign file was generated using the following command.
Peak calling. Peaks were identified using MACS2 (version 2.2.6) with the following parameterization 36 . The input sample was used as the control for the BRD4 ChIP data; the H3K27ac data was processed without an input with MACS2 determining the control by default processes. The BRD4 data was processed with the following parameters -g hs, -p 1e-2,-nomodel,-extsize 121, -B. The H3K27ac data was processed with the parameter set -g hs, -p 1e-2,-nomodel,-extsize 218, -B.
Determination of the final peak set. The called peaks were then intersected with all genes in the hg19 human genome, using bedtools intersect (1 bp overlap) and overlapping regions were removed (https://bedtools.readthedocs.io/en/latest/). The remaining peaks from H3K27ac and BRD4 that did not overlap genes were then intersected using bedtools (1 bp overlap) and regions with both an H3K27ac and BRD4 peak were kept (using the BRD4 coordinates).
Creating BigWigs. The fold enrichment of the bam files was calculated across these peaks for both H3K27ac and BRD4 using macs2 bdgcmp and the -m ppois parameter. As the H3K27ac had no input we felt that in order to allow for fair comparison both H3K27ac and BRD4 bedgraph files should use the -m ppois parameter (we did also generate a fold enrichment aka FE bedgraph for BRD4 to ensure it was comparable to the ppois version). These bedgraph files were then converted to bigwigs using bedGraphToBigWig (version v4) from UCSC.
Calling super enhancers. Super enhancers were then identified using the ROSE 41 pipeline (version 0.1, python 2.7) with default parameters.
Meta-analysis. Meta plots and heatmaps for these data were created using Deeptools (version 3.1.0). We generated matrices using signals from the bigwig files and the overlapping 12,339 peaks as the regions. These matrices were then used for plotting.
ChIP-Seq (H3K9me3). ChIP-Seq analysis for H3K9me3 was performed by Active Motif following their spike in protocol. The following is a modified excerpt from the workflow provided to us.
Sequence analysis. The 75-nt single-end (SE75) sequence reads generated by Illumina sequencing (using NextSeq 500) were mapped to the genome using the BWA algorithm ("bwa aln/samse" with default settings). Alignment information for each read is stored in the BAM format. Only reads that pass Illumina's purity filter, align with no more than 2 mismatches, and map uniquely to the genome were used in the subsequent analysis. In addition, duplicate reads ("PCR duplicates") were removed.
Determination of fragment density. Since the 5´-ends of the aligned reads (="tags") represent the end of ChIP/IP-fragments, the tags were extended in silico (using Active Motif software) at their 3´-ends to a length of 200 bp, which corresponds to the average fragment length in the size-selected library. To identify the density of fragments (extended tags) along the genome, the genome was divided into 32-nt bins and the number of fragments in each bin is determined. This information ("signal map"; histogram of fragment densities) is stored in a bigWig file. bigWig files also provide the peak metrics in the Active Motif analysis program described below.
Peak finding. The generic term "Interval" is used to describe genomic regions with local enrichments in tag numbers. Intervals are defined by the chromosome number and a start and end coordinate. The peak caller used at Active Motif for this project was SICER 91 . This method was used to detect significant enrichments in the ChIP/IP data file when compared to the Input data file or relative to neighboring background regions.
Additional Analysis Steps:. (a) Standard Normalization: In the default analysis, the tag number of all samples (within a comparison group) is reduced by random sampling to the number of tags present in the smallest sample.
(b) Spike-in Adjustment: Spike-in of Drosophila chromatin was performed; the number of test tags was adjusted (again by down-sampling) by a factor that would result in the same number of spike-in Drosophila tags for each sample.
Merged region analysis. To compare peak metrics between 2 or more samples, overlapping Intervals (orange bars in diagram below) were grouped into "Merged Regions" (green bars), which are defined by the start coordinate of the most upstream Interval and the end coordinate of the most downstream Interval (=union of overlapping Intervals; "merged peaks"). In locations where only one sample has an Interval, this Interval defines the Merged Region. The use of Merged Regions was necessary because the locations and lengths of Intervals are rarely exactly the same when comparing different samples. Furthermore, with this approach fragment density values could be obtained even for samples for which no peak was called.
Annotations. After defining the Intervals and Merged Regions, their genomic locations along with their proximities to gene annotations and other genomic features are determined. In addition, average and peak (i.e. at "summit") fragment densities within Intervals and Merged Regions were compiled.
Differential binding analysis. DESeq2 was used to determine regions of differential binding.

Hi-C
In situ Hi-C. OVCAR3 cells were grown under recommended culture conditions in RPMI media supplemented with 10% FBS and 1% penicillin/streptomycin. Four to five million cells were fixed with 1% formaldehyde for 10 min, then cell pellets were flash frozen and stored at −80°C.
In situ Hi-C was performed as previously described 92 . In short, pellets were lysed in ice-cold Hi-C lysis buffer (10 mM Tris-HCl pH 8.0, 10 mM NaCl, 0.2% IGEPAL CA630) with 50 μL of protease inhibitors for 15 min on ice. Cells were pelleted and washed once more using the same buffer. Pellets were resuspended in 50 μL of 0.5% SDS and incubated at 62°C for 7 min. Reactions were quenched with 145 μL water and 25 μL 10% Triton X-100 at 37°C for 15 min. Chromatin was digested overnight with 25 μL of 10X NEBuffer2 and 100U of MboI at 37°C with rotation.
Samples were cooled to RT and 1.6× volumes of pure ethanol and 0.1× volumes of 3 M sodium acetate pH 5.2 were added to each sample, which were subsequently incubated at −80°C for over 4-6 h. Samples were spun at max speed at 2°C for 15 min and washed twice with 70% ethanol. The resulting pellet was dissolved in 130 μL of 10 mM Tris-HCl pH 8.0 and incubated at 37°C for 1-2 h. Samples were stored at 4°C overnight.
DNA was sheared using the Covaris LE220 (Covaris, Woburn, MA) to a fragment size of 300-500 bp in a Covaris microTUBE. DNA was transferred to a fresh tube and the Covaris microTUBE was rinsed with 70 μL of water and added to the sample. A 1:5 dilution of DNA was run on a 2% agarose gel to verify successful shearing.
Sheared DNA was size selected using AMPure XP beads. 0.55× volumes of 2× concentrated AMPure XP beads were added to each reaction and incubated at RT for 5 min. Beads were reclaimed on a magnet and the supernatant was transferred to a fresh tube. 30 μL of 2× concentrated AMPure XP beads were added and incubated for 5 min at RT. Beads were reclaimed on a magnet and washed with fresh 70% ethanol. Beads were dried for 5 min at RT prior to DNA elution in 300 μL of 10 mM Tris-HCl pH 8. Undiluted DNA was run on a 2% agarose gel to verify successful size selection between 300 and 500 bp.
150 μL of 10 mg/mL Dynabeads MyOne Streptavidin T1 beads were washed with 400 μL of 1× Tween washing buffer (TWB; 250 μL Tris-HCl pH 7.5, 50 μL 0.5 M EDTA, 10 mL 5 M NaCl, 25 μL Tween 20, 39.675 μL water). Beads were then resuspended in 300 μL of 2X binding buffer (500 μL Tris-HCl (pH 7.5), 100 μL 0.5 M EDTA, 20 mL 5 M NaCl, 29.4 mL water), added to the DNA sample, and incubated at RT for 15 min with rotation. DNA-bound beads were then washed twice with 600 μL of 1X TWB at 55°C for 2 min with shaking. Beads were resuspended in 100 μL 1× NEBuffer T4 DNA ligase buffer, transferred to a new tube, and reclaimed. Sheared ends were repaired by resuspending the beads in 88 μL of 1× NEB T4 DNA Ligase Buffer with 1 mM ATP, 2 μL of 25 mM dNTP mix, 5 μL of 10 U/μL NEB T4 PNK, 4 μL of 3 U/μL NEB T4 DNA polymerase I, and 1 μL of 5 U/μL NEB DNA polymerase 1, large (Klenow) fragment and incubating at RT for 30 min. Beads were washed two more times with 1× TWB for 2 min at 55°C with shaking. Beads were washed once with 100 μL of 1× NEBuffer 2, transferred to a new tube, and resuspended in 90 μL of 1X NEBuffer 2, 5 μL of 10 mM dATP, and 5 μL of NEB Klenow exo minus, and incubated at 37°C for 30 min. Beads were washed two more times with 1× TWB for 2 min at 55°C with shaking. Beads were washed in 100 μL 1× Quick Ligation Reaction Buffer, transferred to a new tube, reclaimed, and resuspended in 50 μL of 1× NEB Quick Ligation Reaction Buffer. 2 μL of NEB DNA Quick Ligase and 3 μL of an appropriate Illumina indexed adapter (TruSeq nano) were added to each sample before incubating at RT for 15 min. Beads were reclaimed and washed twice with 1× TWB for 2 min at 55°C. Beads were washed in 100 μL 10 mM Tris-HCl pH 8, transferred to a new tube, reclaimed and resuspended in 50 μL of 10 mM Tris-HCl pH 8.
Hi-C libraries were amplified directly off T1 beads with 10 cycles in 5 μL of PCR primer cocktail, 20 μL of Enhanced PCR mix, and 25 μL of DNA on beads. The PCR settings were as follows: 3 min at 95°C followed by 4-12 cycles of 20 s 98°C, 15 s at 60°C, and 30 s at 72°C. Samples were held at 72°C for 5 min before lowering for holding at 4°C. Amplified samples were transferred to a new tube and brought to 250 μL in 10 mM Tris-HCl pH 8.
Beads were reclaimed and the supernatant containing the amplified library was transferred to a new tube. Beads were resuspended in 25 μL of 10 mM Tris-HCl pH 8 and stored at −20°C. 0.7× volumes of warmed AMPure XP beads were added to the supernatant sample and incubated at RT for 5 min. Beads were reclaimed and washed once with 70% ethanol without mixing. Ethanol was aspirated. Beads were resuspended in 100 μL of 10 mM Tris-HCl pH 8, 70 μL of fresh AMPure XP beads were added, and the solution was incubated for 5 min at RT. Beads were reclaimed and washed twice with 70% ethanol without mixing. Beads were left to dry and DNA was eluted in 25 μL of 10 mM Tris-HCl pH 8. The resulting libraries were next quantified by Qubit and Tapestation. A low-depth sequence was performed first using the Miniseq sequencer system (Illumina) and analyzed using the Juicer pipeline to assess quality. The resulting libraries underwent paired-end 2 × 150 bp sequencing on an Illumina NovaSeq sequencer. Each replicate was sequenced to an approximate depth of 730 million reads. The full sequencing depth was approximately 2.92 billion reads.
Hi-C data processing and analysis. In situ Hi-C datasets were processed using dietJuicer, a modified version of the Juicer Hi-C pipeline (https://github.com/ EricSDavis/dietJuicer), using default parameterization 93 . Reads were aligned to the hg19 human genome (using Mbol restriction enzyme) with bwa (version 0.7.17). Four biological replicates were aligned and then these replicates were merged for a total of 2,922,558,308 Hi-C read pairs in OVCAR3 cells yielding 2,598,024,810 valid Hi-C contacts (88.90%). The resulting Hi-C contact matrix was next normalized with the "KR" matrix balancing algorithm. This was done in order to adjust for regional background differences in chromatin accessibility and allow for proper visualization of this data 94 .
CRISPR-KO gene targets were identified as direct or indirect targets using Hi-C contact frequency. Specifically, we compared the fold-change in observed over expected contact frequency between SE14 or SE60 and their respective gene targets with 100 permutations of distance-matched region-gene pairs as controls. Since distance-matching is only relevant for regions within a chromosome, we restricted our analysis to intra-chromosomal pairs. Direct targets were defined as SE-gene pairs with an observed/expected contact frequency greater than the 75th percentile of the control distribution. We performed this analysis on (1)  Hi-C ABC analysis. We used the ABC method 62 with slight modifications to determine direct super-enhancer (SE) KO target genes. BRD4 signal for each SE was used to measure enhancer activity and 50 kb resolution, KR-normalized Hi-C counts were used to measure contact between each SE and their putative target genes. Putative targets were defined as expressed genes (baseMean > 100 counts). Activity and contact vectors were log2-transformed and scaled before being multiplied together to create an ABC score for each SE-target pair.
Single-cell analysis Data acquisition. We obtained the single-cell RNA-seq and single-cell ATAC-seq data from the GEO accession number GSE173682.
scRNA-seq data processing and barcode quality-control (QC). The filtered featurebarcode matrix was converted into a Seurat object for each patient tumor sample using the Seurat R package (Seurat version 3.2) 97 . QC and doublet removal were performed for each patient dataset individually to emrich for high-quality cells. Outlier cells were defined in each of the following metrics: log(UMI counts) (>2 MADs, low end), log(percent mitochondrial read count +1) (>2 MADs, high end), and log(number of genes expressed) (>2 MADs, low end). Non-outlier cells, according to all three criteria, were kept for doublet detection. Cells marked as doublets by both DoubletDecon 98 (version 1.1.5) and DoubletFinder 99 (version 2.0.3) were removed from downstream analysis.
scRNA-seq clustering and cell type annotation. Seurat objects were normalized using Seurat's NormalizeData() with the normalization method set to "LogNormalize." Seurat's FindVariableFeatures() was used for feature selection with the selection method set to "vst" and the number of top variable features set to 2000. Prior to principal component analysis (PCA), gene expression values were scaled using Seurat's ScaleData(). The top 2000 most variably expressed genes were summarized by PCA and the cells were visualized in a two-dimensional UMAP plot using Seurat's RunUMAP() with 50 principal components (PCs), as suggested by the results of Seurat's JackStraw(). To cluster cells, graph-based Louvain clustering was performed using Seurat's FindNeighbors() with all 50 PCs and Seurat's FindClusters() with a resolution of 0.7. scRNA-seq UMAP plots were generated in R using ggplot2.
Cell type annotation was performed using the R package SingleR 100 and was verified with gene signature enrichment scores using Seurat's AddModuleScore(). Both scRNA-seq datasets used in this study were annotated based on a reference scRNA-seq dataset from a human ovarian tumor (sample ID: HTAPP-624-SMP-3212) 101 . The individual patient datasets were then combined using Seurat's merge() and subsequently reprocessed according to the normalization, feature selection, and clustering methods described above. The resulting clusters in the merged dataset were annotated based on the majority of cell type label within each cluster. SingleR cell type annotations were verified by calculating cell type gene signature enrichment scores from PanglaoDB 102 using Seurat's AddModuleScore(). scRNA-seq differential gene expression analysis. Differential gene expression was computed using Seurat's FindMarkers() with the "test.use" parameter set to "wilcox" for the Wilcoxon Rank Sum test. Genes with a Bonferroni-corrected Pvalue ≤ 0.01 and average logFC ≥ 0.1 were deemed upregulated in the cancer epithelial fraction relative to the remaining cell type clusters.
scATAC-seq data processing and barcode quality-control (QC). The scATAC-seq fragments file for each patient tumor sample was read into the R package ArchR (version 0.9.3) to perform barcode quality control and doublet removal 103 . To enrich for cellular barcodes, log10(TSS enrichement+1) and log10(number of unique fragments) were used as QC metrics, both of which showed a bimodal distribution. Gaussian mixture models (GMM), implemented in the R package mclust 104 , were used to estimate barcode cutoff thresholds for log10(TSS enrich-ment+1) and log10(number of unique fragments). Barcodes above these GMMestimated thresholds in both metrics were retained. ArchR's addDoubletScores() function was used to calculate doublet enrichment scores with the knnMethod set to "UMAP" and putative doublets were then filtered out using ArchR's filter-Doublets() with default parameters.
scRNA-seq cell type label transfer to scATAC-seq. Prior to transferring labels from scRNA-seq to scATAC-seq, ArchR's addGeneScoreMatrix() was used to infer gene activity scores in scATAC-seq using default parameters. ArchR's addGeneInte-grationMatrix() was used to assign each of the scATAC-seq cells a cell type subcluster identity from the matching scRNA-seq data and an associated label prediction score 97 . Of note, this label transfer procedure was constrained to only align cells of the same patient dataset. Only scATAC-seq cells with a label prediction score >0.5 were included in downstream analyses. Moreover, only inferred cell type subclusters with >30 cells were included in the downstream analysis.
scATAC-seq peak calling and data visualization. Pseudo-bulk replicates were generated for each inferred cell type subcluster using the R package ArchR and pseudo-bulk peak calling was performed within each inferred cell type subcluster using MACS2 36,103 . The peak calls from each inferred cell type subcluster were then merged into a single profile using ArchR's default iterative overlap procedure to create a merged peak by barcode matrix across all cellular barcodes from both patient tumor samples. ArchR's plotBrowserTrack() was used to visualize the scATAC-seq coverage per inferred cell type. scATAC-seq differential peak accessibility for determining cancer-enriched enhancers. ArchR's getMarkerFeatures() was used for differential peak accessibility analysis with the bias argument set to include both "TSSEnrichment" and "log10(number of fragments)." Differentially accessible peaks (DEPs) were identified for each cell cluster by comparing the accessibility values of peaks across all cells in a cluster (group 1) relative to the accessibility values for a group of background cells matched for TSS enrichment and read depth (group 2). This comparison was made between cancer clusters (group 1) and all remaining cell-type clusters (group 2). Peaks with Benjamini-Hochberg FDR ≤ 0.10 and Log2FC ≥ 0.25 were deemed cancer-enriched with statistically significant increased accessibility relative to the non-cancer fraction.
Enhancer motif analysis in scATAC-seq. Bedtools getfasta() was used to extract the sequences of select cancer-enriched enhancers according to the hg38 reference genome 105 . FIMO motif scanning with default parameters was applied to the enhancer sequences using a motif database supplied by JASPAR2020 64,106 . FIMO motif results were ranked by Benjamini-Hochberg corrected q-values and TF expression in the cancer fraction by summing the normalized TF counts across all cells within the cancer epithelial clusters. Seurat's VlnPlot() was used to generate TF expression violin plots.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
Data generated in this study (CRISPRi/CRISPR KO RNA-seq, Hi-C, and H3K9me3 ChIP-Seq) have been uploaded and are publicly available in the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/) database under the accession number GSE174259.
The single cell genomics data were downloaded from GEO accession number GSE173682 43 . The H3K27ac ChIP-seq and BRD4 ChIP-seq were downloaded from GEO accessions GSE101408 28 and GSE77568 28 , respectively. All FTSEC data was downloaded from GEO accession GSE68104 38 . TCGA expression and copy-number data 10 (RNA-seq and CNV data) was downloaded using the TCGA repository Firebrowse (http:// firebrowse.org/) and survival data was attained from the supplement of the TCGA clinical paper 44 . All ENCODE data was downloaded from the Screen database 37 (https:// screen.encodeproject.org/#).
Kaplan-Meier plots for gene set survival analysis were created with the publicly available KM Plot tool 34,35 (https://kmplot.com/analysis/). Patient data from Fig. 1b, c are publicly available through CBioPortal 6,31 (https://bit.ly/3QY91sa). The bar chart from Fig. 1b can be found under Cancer Types Summary. The boxplot from Fig. 1c can be generated under Plots by plotting TCGA PanCanAtlas Cancer Type Acronym vs. mRNA Expression, RSEM (Batch normalized from Illumina HiSeq_RNASeqV2) (log2(value + 1)) and sorting the categories by a median. For both Fig. 1b, c, the 16 highest altered/expressed TCGA cancer types are presented.
The remaining data are available within the article, Supplementary Information, or Source Data file. Source data are provided with this paper.